Approche variationnelle pour le calcul 
bayésien dans les problèmes inverses en 

imagerie 

Ali Mohammad-Djafari 

'■ 

O ; Laboratoire des signaux et systèmes 

^ : (UMR 08506, CNRS-SUPELEC-Univ Paris Sud) 

^ ; Supélec, Plateau de Moulon, 91192 Gif-sur- Yvette Cedex, France 

l> 



> 

00 



o 



27 avril 2009 



Résumé 



Dans une approche bayésienne non supervisée pour la résolution d'un problème 
inverse, on cherche à estimer conjointement la grandeur d'intérêt / et les para- 
mètres 6 à partir des données observées g et un modèle A4 liant ces grandeurs. Ceci 
se fait en utilisant la loi a posteriori conjointe p{f ,6\g; M). L'expression de cette 
loi est souvent complexe et son exploration et le calcul des estimateurs bayésiens 



^ ■ nécessitent soit les outils d'optimisation de critères ou de calcul d'espérances des 

On . lois multivariées. Dans tous ces cas, il y a souvent besoin de faire des approxima- 



tions. L'approximation de Laplace et les méthodes d'échantillonnage MCMC sont 



^ . deux approches classiques (analytique et numérique) qui ont été explorés avec suc- 

' cès pour ce fin. Ici, nous étudions l'approximation de p{f, 0\g) par une loi séparable 

en / et en 0. Ceci permet de proposer des algorithmes itératifs plus abordables 
en coût de calcul, surtout, lorsqu'on choisit ces lois approchantes dans des familles 
des lois exponentielles conjuguées. Le principal objet de ce papier est de présenter 
les différents algorithmes que l'on obtient pour différents choix de ces familles. À 
titre d'illustration, nous considérons le cas de la restauration d'image par déconvo- 
lution simple ou myope avec des a priori séparables, markoviens simples ou avec 
des champs cachés. 
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Abstract 

In a non supervised Bayesian estimation approach for inverse problems in imaging 
Systems, one tries to estimate jointly the unknown image pixels / and the hyperparameters 
given the observed data g and a modcl M. linking thcsc quantitics. This is, in général, 
done through the joint posterior law p{f ,d\g] M). The expression of this joint law is 
often very complex and its exploration through sampling and computation of the point 
estimators such as MAP and posterior means need either optimization of or intégration of 
multivariate probability laws. In any of thèse cases, we need to do approximations. Laplace 
approximation and samphng by MCMC are two approximation methods, respectively 
analytical and numerical, which have been used before with success for this task. In 
this paper, we explore the possibility of approximating this joint law by a separable one 
in / and in d. This gives the possibility of developing itérative algorithms with more 
reasonable computational cost, in particular, if the approximating laws are choosed in 
the exponential conjugate families. The main objective of this paper is to give détails of 
différent algorithms we obtain with différent choices of thèse families. To illustrate more 
in détail this approach, we consider the case of image restoration by simple or myopie 
deconvolution with separable, simple markovian or hidden markovian models. 
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1 Introduction 



Une présentation simplifiée et synthétique des problèmes inverses en imagerie, en se 
plaçant en dimensions finies, consiste à vouloir retrouver une grandeur inconnue / à partir 
des observations g d'une grandeur observée, en supposant connaître un modèle M. qui les 
lient. La forme la plus simple de ce modèle est un modèle linéaire de la forme 

M: g = Hf + e (1) 

où on suppose que toutes les erreurs de modélisation et de mesure peuvent être représentées 
par e. Notons aussi que g et f sont, en général, des vecteurs de grandes dimensions, 
ce qui signifie que nous considérons ici le cas discrétisé où g contient l'ensemble des 
grandeurs mesurées et / l'ensemble des valeurs qui décrivent la grandeur inconnue. Dans 
ce contexte H est une matrice dont les éléments sont définies par le modèle et les étapes 
de discrétisation du problème. 

Dans une approche estimation bayésienne non supervisée pour résoudre un problème 
inverse, d'abord on utihse ce modèle pour définir la loi de probabilité p{g\f,Oi;M.) où 
6i représente l'ensemble des paramètres qui décrivent cette loi. Lors que cette fonction 
est considérée comme une fonction de / et de 6i, elle est appelée la vraisemblance des 
inconnues / et du modèle Ai. Son expression s'obtient à partir de la loi de probabilité 
des erreurs e en utilisant le modèle ([T]). Par exemple, lorsque e dans ce modèle est modélisé 
par un vecteur aléatoire centré, blanc, gaussien et de covariance fixée = aj^l, on a 

pig\f,e,;M)=MiHf,^,) (2) 

où Oi = (j^ . D'autres lois avec d'autres paramètres d\ peuvent bien sûr être utilisées. 

La deuxième étape dans cette approche est l'attribution ou le choix d'une loi dite a 
•priori p{f\02; Ai) pour les inconnues /, où 62 représente ses paramètres. La troisième 
étape consiste à écrire l'expression de la loi a posteriori des inconnues / : 

,,,„„,,, p{g,mM) p{g\f,er,M)p(m;M) 

où on suppose implicitement connaître l'ensemble des paramètres 6 = {61, 62). Mais, dans 
un cas réel, nous sommes amenés souvent à les estimer aussi. Pour cela, dans l'approche 
bayésienne, on leur attribue aussi une loi a priori p{6\Ai), et l'on obtient alors une loi a 



posteriori conjointe des inconnues / et des hyperparamètres 6 = {61,62) 

p{gJ,6\M) 



p{f,6\g;M) = , , . 

p{g\M) 

p{g\f,6,;M)p{m;M)p{e\M) 



P{9\M) 

Dans cette relation, le dénominateur 

P(9\M) = Il p(g\f, 6- M) p{f\6- M) p{e\M) df de (5) 

est la vraisemblance marginale du modèle M. dont son logarithme lnp{g\M.) est appellé 
évidence du modèle Ai. 

Afin d'introduire les notions qui vont être utilisées dans la suite de ce travail, il est 
intéressant de mentionner que, pour n'importe quelle loi de probabilité q{f , 6) (dont nous 
verrons le choix et l'utilité par la suite), l'évidence du modèle lnp(p|A1) vérifie 



\M9\M) = \nJJp{g,f,e\M)dfde 



(d'après l'inégalité de Jensen : ln(E{p/g}) > E{ln(p/g)}). Aussi, notant par 

TM = jj,(f.0)^'-^^§^éfie (T) 

et par 

KLta : p) = h. -j^(|^d/ de (8) 

on montre facilement (en remplaçant p{g, /, 6\A4) — p{f, e\g; M) p{g\M.) dans l'expres- 
sion de ^(ç)) que 

\iip{g\M) = T{q) + Kh{q:p). (9) 

Ainsi T{q), appelée l'énergie hbre de g(/, 6) par rapport à p{g, /, e\M.), est une hmite 
inférieure de \n.p{g\M.) car KL(ç : p) > 0. Par la suite, nous allons écrire l'expression de 

T{q) par 

T{q) = {\Tip{gJ,6-M)), + n{q) (10) 

011 nous avons utilisé la notation < . >q pour l'espérance suivant la loi q et Ti-iq) est 
l'entropie de q : 

H{q) = - ff q{f,6)\nq{f,6) df d6 (11) 



Arrivé à ce stade, les questions posées sont : 

- Inférence : Etant donnée les expressions des lois a posteriori ([3]) et (jl]), comment 
les utiliser pour définir une solution au problème inverse décrit en ([T]) ? 

- Sélection de modèle : Comment peut-on sélectionner un modèle parmi un en- 
semble de modèles A^j. 

En ce qui concerne le problème de l'inférence, les principaux choix sont les estimateurs 
au sens du Maximum a posteriori (MAP) ou au sens de la moyenne a posteriori (PM). 
Dans le premier cas, on a besoin des outils d'optimisation et dans le deuxième cas des 
outils d'intégration (analytique ou numérique). Pour la sélection du modèle, nous nous 
contenterons ici de noter que l'expression de la vraisemblance du moèle ([5]) peut être 
utilisés à cette fin. Dans ce papier, nous nous focalisons sur la première question où on 
cherche à inférer f et utilisant la loi a posteriori jointe (jl]). 

Pour les problèmes inverses, la solution au sens du MAP a été utilisée avec succès 
pour sa simplicité et en raison de son lien avec l'approche déterministe de la régulari- 
sation. Mais, il y a des situations où cette solution ne donne pas satisfaction, et où la 
solution au sens de la moyenne a posteriori peut être préférée. Mais, les situations où on 
peut avoir une solution analytique pour les intégrations qui sont nécessaires pour obtenir 
ces estimées sont rares. Il y a alors pratiquement deux voies : 

Intégration numérique par échantillonnage : Il s'agit d'approcher les espérances 
par des moyennes empiriques des échantillons générés suivant la loi a posteriori . Toute 
la difficulté est alors de générer ces échantillons, et c'est là qu'interviennent les méthodes 
de MCMC (Markov Chain Monté Carlo). Le principal intérêt de ces méthodes est qu'elles 
permettent d'explorer l'ensemble de l'espace de la loi a posteriori , mais l'inconvénient 
majeur est leur coût de calcul qui est dû au nombre important d'itérations nécessaire 
pour la convergence des chaînes et le nombre important de points qu'il faut générer pour 
obtenir des estimations de bonnes qualités. 

Approximation de la loi a posteriori par des lois plus simples : Il s'agit de 
reporter le calcul des intégrales après une simplification par approximation de la loi a 
posteriori . Une première approximation utilisée historiquement est Approximation de 
Laplace où on approxime la lois a posteriori par une loi gaussienne. Dans ce cas, les deux 
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estimateurs MAP et PM sont équivalents et tous les calculs deviennent analytiques. 

Une deuxième solution est d'approcher la a posteriori par une loi séparable, ce qui 
permet de réduire la dimension des intégrations. Cette voie est plus récente et le principal 
objet de ce papier. 

De façon générale, l'idée d'approcher une loi conjointe p{x) de plusieurs variables x 
par une loi séparable q{x) = Yljlji^j) n'est pas nouveau et peut être trouvée dans la 
littératures de la fin des années 90 : [UElIllIlllHlEllT]. Le choix d'un critère pour me- 
surer la qualité de cette approximation et l'étude des effets de cette approximation sur 
les qualités des estimateurs obtenus apparaît dans les travaux plus récents [H [9l [101 [Hl 
[T2| [T3l [m [ini [m [TTI [m [H]. L'usage de cette approche en estimation des paramètres 
d'un modèle d'observation avec des variables cachées en statistique est également ré- 
cente [2ni[ï9l[2ll[23[23l[2a[25l[26l[27l[2a Dans la plupart de 
ces travaux, l'application est plutôt en classification en utilisant un modèle de mélanges. 
Dans le domaine du traitement du signal, ils utilisent des modèles de mélange avec des 
étiquettes des classes modélisées par des chaînes de Markov. Dans le domaine du trai- 
tement d'image, la plupart de ces travaux sont consacrés à la segmentation d'image en 
utilisant soit un modèle séparable ou markovien pour les étiquettes. Cette approche a été 
aussi utilisé récemment en séparation de sources [361 EU EH] et en traitement des images 
hyperspectrales [391110]. L'application réelle de cette approche pour simplifier les calculs 
bayésiens dans les problèmes inverses avec un opérateur mélangeant est l'originalité de ce 
travail [Hl |12|, |13]. En effet, dans la plupart de ces travaux, avec les notations utilisées 
dans ce papier, on suppose qu'on a observé directement / dont la loi est modéhsée par 
un mélange de gaussiennes et le principal objectif de ces méthodes est l'estimation des 
paramètres de ce mélange et la sélection du modèle. Dans d'autres travaux on considère 
un modèle d'observation ponctuel où Çi = h{fi) + et l'objective est la segmentation de 
l'image /. Les travaux dans lesquels, on utilise l'approche variationnelle pour des pro- 
blèmes inverses de restauration ou de reconstruction d'images sont assez rare. On trouve 
essentiellement les travaux [Hl HÏÏl |16] qui considèrent le cas des problèmes inverses li- 
néaires, mais soit avec des modèles a priori de Gauss-Markov ou avec des variables cachées 
de contours ou d'étiquettes des régions séparables. La dépendance spatiale de ces variables 
cachées n'est pas prise en compte. 
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La suite de ce papier est organisée de la manière suivante : une présentation très syn- 
thétique de l'approche variationnelle est fournie dans la section 2. Dans la section 3, nous 
nous intéresserons à l'approximation de la loi a posteriori p{f, 6\g; A4) pour les problèmes 
inverses. Il s'agit alors d'appliquer l'approche variationnelle pour ce cas particulier avec 
différent choix pour les familles des lois approchantes. Dans la section 4, nous détaillerons 
l'application de la méthode au cas de la restauration simple ou myope d'image. Enfin, dans 
la section 5, nous décrirons la manière d'appliquer cette méthode au cas de la restauration 
d'image avec une modélisation a priori plus complexe : un modèle hiérarchique de Gauss- 
Markov avec un champ caché de Potts pour des étiquettes des régions dans l'image. En 
effet, cette modélisation convient pour bien des cas de restauration et de reconstruction 
d'images dans des applications oii l'image recherchée représente un objet composé d'un 
nombre fini de matériaux homogènes. C'est exactement cette information a priori qui est 
modélisée par un modèle de mélange avec des étiquettes z markoviennes [17]. Ici, donc, 
la loi a posteriori p{f,z,6\g;Ai) doit être approximée par des lois plus simples à utili- 
ser. Finalement, dans la section 6, nous résumons les apports de ce papier. Il est à noter 
cependant, que dans ce papier, seules les principes des méthodes proposées sont présentés 
et les détails de mise en oeuvre de ces méthodes, ainsi que des résultats de simulations et 
évaluation des performances de ces algorithmes sont reportées à un prochaine papier. 

2 Principe de l'approche variationnelle 

Considérons le problème général de l'approximation d'une loi conjointe p{x\Ai) de 
plusieurs variables x par une loi séparable q{x) = YljÇij{xj). videment, il faut choisir un 
critère. Considérons le critère de divergence KL(g : p) entre p et q : 



et cherchons la solution qui le minimise. La solution optimale ne peut être calculée 
que d'une manière itérative, car ce critère n'étant pas quadratique en g, sa dérivée par 
rapport à q n'est pas linéaire et une solution explicite n'est pas disponible. On peut alors 
envisager une solution itérative coordonnée par coordonnée. 

Notons q^j{x_j) = Yli^éj qi{xi), ce qui permet d'écrire q{x) = qj{xj) q^j{x_j) où x_j = 
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{xi,i 7^ j} et développons ce critère : 

KL(g:p) = - f q(x)\n^^^$^ dx 
= -Hiq)-{lnpix\M))^^^^ 

j 

= -Y^H{q,)~ I {\r.p{x\M))^_^q,{x^) àx, (12) 

où H{q) est l'entropie de q et H{qj) est l'entropie de qj. 

Notons que si q-j est fixée KL (g : p) est convexe en qj et sa minimisation sous la 
contrainte de normalisation de qj s'écrit 



{\np{x\M))^_^ (13) 



avec 



Cj = / exp 



{\np{x\M))g_ dxj (14) 



On note que l'expression de qj dépend de celle de q-j et que, l'obtention de qj (étant 
donnée g_j) se fait aussi d'une manière itérative en deux étapes : 



Q(g,) = H{q,) + {\npix\M))^^^ 



^f^^^ = arg maxg^ { Q{qj)} , 
ce qui ressemble à un algorithme du type EM (Espérance-Maximisation) généralisé, car 
dans la première équation nous avons à calculer une espérance et dans la deuxième étape 
nous avons une maximisation. 

Les calculs non paramétriques sont souvent trop coûteux. On choisit alors une forme 
paramétrique pour ces lois de telle sorte que l'on puisse, à chaque itération, remettre à jour 
seulement les paramètres de ces lois, à condition cependant que ces formes ne changent pas 
au cours des itérations. La famille des lois exponentielles conjuguées ont cette propriété 

La première étape de simplification de ce calcul est donc de considérer la famille des 
lois paramétriques qj{xj\6j) oii Oj est un vecteur de paramètres. En effet, dans ce cas, 
l'algorithme itératif précédent se transforme à : 

Qid,) =H{q,ix,\e,)) + {\npix\M))^^^,^^,^), ^^^^ 
ef+'^ =argmax«^,{Q(0,)}, 



(15) 



ce qui ressemble plus à un algorithme du type EM. 

Une deuxième étape de simplification est de choisir pour p{x\0; A4) la famille des lois 
exponentielles conjuguées : 

p{x\e; M) = f{x) g{e) exp [(t>{efu{x)] (17) 

où <p{0) est le vecteur de paramètres naturel et f{x), g{6) et u{x) sont des fonctions 
connues. Il est alors facile de montrer que qj{xj) dans l'équation f|T3l) sera aussi dans la 
famille des lois exponentielles conjuguées, et par conséquence q{x\0) restera dans la même 
famille et nous aurons juste à remettre à jours ses paramètres. 

Remarque : Cette famille de lois a une propriété dite conjuguée dans le sens que si on 
choisit comme a priori 

p{e\r], u) = h{r], u) giey exp [iy^<t>{e)] (18) 

l'a posteriori correspondant 

p{0\x; M) oc p{x\0; M) pi0\ri,u) (19) 

sera aussi dans la même famille (ÎT8|) . La famille des lois a priori (fT7|) est alors dites 
conjuguée de la famille des lois (ÎT6|) . 

Pour cette famille de lois, on a 

{\np{x\9; M))^ = {\nf{x))^ + \ng{e) + 0(0)^ {u{x)}^ (20) 

et donc la forme de la loi séparable qj correspondante est 

q,{x,\e) oc g{è)exp [< f{x) +^{ef < u{x) (21) 

où désigne les paramètres particuliers de qj qu'il faut mettre à jour au cours des itéra- 
tions. 

Par ailleurs, sachant que q{x) est séparable, si lnp{x\0; M) est polynomial en x, le 
calcul de l'espérance (12U|) . mais surtout celui de sa dérivée par rapport aux qj et aux 
paramètres 0, qui sont nécessaire pour l'optimisation, seront facilités. Comme nous le 
verrons par la suite, pour les problèmes inverses, nous avons choisi ce genre de lois. 
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3 Approche variât ionnelle pour les problèmes inverses 
linéaires 



Nous allons maintenant utiliser ces relations pour décrire le principe de l'approche 
variationnelle au cas particulier des problèmes inverses ([î]) où l'utilisation directe de la 
loi a posteriori conjointe 6\g., M) dans l'équation (jl]) est souvent trop coûteuse pour 
pouvoir être explorée par échantillonnage direct de type Monté Carlo ou pour calculer les 
moyennes a posteriori 

f = JJ fp{f,0\g;M)dedf (22) 

et 

= [ [ 9p{f,0\g;M)àfdO. (23) 



En effet, rare sont les cas où on puisse trouver d'expressions analytiques pour ces inté- 
grales. De même l'exploration de cette loi par des méthodes de Monté Carlo est aussi 
coûteuses. On cherche alors à l'approcher par une loi plus simple q{f, 6). Par simplicité, 
nous entendons par exemple une loi q qui soit séparable en / et en : 

q{f,0)=q^{f)q2{9) (24) 

videmment, cette approximation doit être faite de telle sorte qu'une mesure de distance 
entre q et p soit minimale. Si, d'une manière naturelle, on choisi KL (g : p) comme cette 



(Îl, $2) = arg min {KL(gig2 ■ p)} = arg max {J^{qiq2)} (25) 

(91:92) (91,92) 

et sachant que KL(gig2 : p) est convexe en qi à q2 fixée et vice versa, on peut obtenir la 
solution d'une manière itérative : 

qi = argmiug, {KL(çiÇ2 : p)} = argmaxg, {J^(gig2)}; 

(26) 

q2 = argmmgj {KL(gig2 ■ p)} = argmaXg^ {J^(gig2)} • 

Utilisant la relation ([7j), il est facile de montrer que les solutions d'optimisation de ces 
étapes sont 



qi{f) ocexp 
^2(0) oc exp 



{lnp{g,f,e-M))^^^f,^ 
{\np{g,f,e;M))^^^f^ 



(27) 
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Une fois cet algorithme convergé vers ql{f) et q^{6), on peut les utiliser d'une manière 
indépendante pour calculer, par exemple les moyennes f* = Jf Qiif) d/ et 

0* = je ^(6>) de. 

Comme nous l'avons déjà indiqué, les calculs non paramétriques sont souvent trop 
coûteux. On choisit alors une forme paramétrique pour ces lois de telle sorte que l'on 
puisse, à chaque itération, remettre à jour seulement les paramètres de ces lois. Nous 



exammons ici, trois cas : 



3.1 Cas dégénéré 

On prend pour qi{f) et q2{e) des formes dégénérées suivantes : 

=^(/-/) (28) 
q2{e\e) =ô{e-e) 

Par conséquence, au cours des itérations, nous n'aurons qu'à remettre à jour / et e. En 
remplaçant qi{f) et Ç2(0) dans les relations ( !27l) on obtient : 

qi{f)^Pi9,f,0;M)^p{f,è\g;M) 
UO) oc pig, /, e, M) oc Pif, e\g; M). 

Il est alors facile de voir que la recherche de / et au cours des itérations devient 
équivalent à : 



/= argmax^ {p{9, f,0;M)^= argmax/ ^p{f,e\g;M)^ 
e= argmaxe ^p{g,f,e]M)^ = argmax^ A^)} 



(30) 



On remarque alors que l'on retrouve un algorithme de type MAP joint ou ICM (Iterated 
Conditional Mode). L'inconvénient majeur ici est qu'avec le choix (l28l) . les incertitudes 
liées à chacune des inconnues / et ne sont pas prise en compte pour l'estimation de 
l'autre inconnue. 

3.2 Cas particulier conduisant à l'algorithme EM 

On prend comme dans le cas précédent une forme dégénérée pour q2{e) = 5{e — e), 
ce qui donne 

qiif) oc p{g, /, è; M) oc p{f,e\g; M) oc p(/|è, g; M) (31) 
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ce qui signifie que qi{f) est une loi dans la même famille que la loi a posteriori p{f\6, g; A4). 
Évidemment, si la forme de cette loi est simple, par exemple une gaussienne, (ce qui est 
le cas dans les situations que nous étudierons) les calculs seront simples. 

A chaque itération, on aurait alors à remettre à jour qui est ensuite utilisé pour 
trouver qi{f\0) = p{f\0,g;J^), elle même utilisée pour calculer 

Qie,0) = {lnp{g,f,0;M))^^^f^~,y (32) 

On remarque facilement l'équivalence avec l'algorithme EM qui se résume à : 

E : Qie, e)= (lnp(g, /, 6- A<)) , ^^^^ 
M: = argmaxe |g(0,ê)| . 

L'inconvénient majeur ici est que l'incertitude liée à 6 n'est pas prise en compte pour 
l'estimation de /. 

3.3 Choix particulier proposé pour les problèmes inverses li- 
néaires 

Il s'agit de choisir, pour gi(/) et q2{0) des lois dans les mêmes familles de lois que 
celles de p{f\g, 6) et de p{0\g, /). En effet, comme nous le verrons plus bas, dans le cas 
des problèmes inverses linéaires ^ avec des choix approprié pour les lois a priori associée 
à la modélisation directe du problème, ces lois a posteriori conditionnelles restes dans les 
mêmes familles, ce qui permet de profiter de la mise à jour facile de ces loi. 

Dans ce travail, dans un premier temps, nous allons considéré le cas des problèmes 
inverses linéaires ^ : g = H f + on H représente la forme discrétisé de la modélisation 
directe du problème et e représente l'ensemble des erreurs de mesure et de modélisation 
avec des hypothèses suivantes : 



(34) 



p{g\H,f,9,;M) = //{H f , {1/9,)!) 
p{f\6f-M) = M{0,{l/ef){D)Di 

Piôf-M) = Q{afo,l3fo), 

où Df est la matrice des différences finies d'ordre un ou deux et 6 = {de = l/a'^,df = 
l/cr|). On obtient alors facilement les expressions dep{g, f\6] A4), p{f\0, g; Ai) et p{6\g, f]A4) 
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qui sont : 

p{g,f\H,9„M) = M{Hf,{l/e,)I)x 

Af{0,{l/ef){D}Df)-^), 
p(f\g,H,ef,M) = mîlf,%), (35) 
p{0e\9;M) = Ç{ae,fy, 
piOf\9;M) = Q{af3f) 
où les expressions de jî^, Sy, (cèe,/9e) et {af,(3f) sont : 

S/ = [< > < ^/ > D)Df]-^ 

= <fe[^*^ + avec A = |g 

jù^ = < > ê/JT*gr, 

«e = aeo + M/2, (36) 

K = /?eo + 1/2 < e*e >, 

S/ = afQ + N/2, 

Pf = /3/o + l/2Tr{£)p/ <//*>}, 



ou 



<e*e>=[5-if/x^]*[^-lîiù^], (37) 

< >= 0;e//?e, < Of >= Oif/(3f 

L'algorithme de mise à jour ici devient : 

- Initialiser = «eo, Â; = /^eo, = «/o, = /3/o 

- Mettre à jour jusqu'à la convergence : 
< > et < ^/ > puis et puis /Ij 

et puis «e, Pe, OÎf, (5f. 

La principale difficulté ici est l'inversion de la matrice [H^H + AD^I?/] qui est de 
dimensions trop grande. Une solution est de choisir pour q{f\g,0) aussi une loi sépa- 
rable q{f\g,0) = Yljlifjlo^^)- L'autre solution, comme nous la verrons dans la section 
suivante, est d'utiliser la structure spécifique de la matrice à inverser pour trouver un 
algorithme d'inversion convenable. 
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4 Application en restauration d'image 

Dans le cas de la restauration d'image où. H a une structure particulière, et où l'opé- 
ration Hf représente une convolution de l'image / avec la réponse impulsionnelle h, la 
partie difficile et coûteuse de ces calculs est celle du calcul de / qui peut se faire à l'aide 
de la Transformée de Fourier rapide [52] . 

De même, l'approche peut être étendue pour le cas de la restauration aveugle ou myope 
où on cherche à la fois d'estimer la réponse impulsionnelle h, l'image / et les hyperpa- 
ramètres 6. Pour établir l'expressions des différentes lois dans ce cas, nous notons que le 
problème direct, suivant que l'on s'intéresse à / (déconvolution) ou k h (identification de 
la réponse impulsionnelle), peux s'écrire 

g(r) = h{r)* f{r) + e{r) = f{r)*h{r) + e{r) 
g = H f + e = F h + e 
où la matrice H est entièrement définie par le vecteur h et la matrice F est entièrement 
définie par le vecteur /. 

Pour permettre d'obtenir une solution bayésienne pour l'étape de l'identification, nous 
devons aussi modéliser h. Une solution est de supposer h = où la matrice $ est 
une matrice telle que $îi> représente la convolution 0(r) * w{r). Ainsi les colonnes de 
$ représentent une base et les éléments du vecteur w représentent les coefficients de la 
décomposition de h sur cette base. On a ainsi 

g(r) = {4>*w)* f{r) + e{r) = f *{4>*w){r) + e{r) 

[d'à) 

g = ^Wf + e = F^w + e 

où la matrice W est entièrement définie par le vecteur w. 

Le problème de la déconvolution aveugle se ramène à l'estimation de / et îi> avec des 

lois 

p{g\w, /, S,) = ATi^Wf, S,) = f^{F^w, S,) (40) 



avec 



= diae 



et p{6ej = Qiaeo,f3eo), 
p{f\ef)=^{0,{ejD}Df)-') (41) 



avec 



p{ôf) = Ç{afo,(3fo) 
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et 



avec 



p{w\ci) = \{p{w^\a^.) = \{N{Q, — ) (42) 

Où-}!)- 



(43) 



3 3 

Avec ces lois a priori , il est alors facile de trouver l'expression de la loi conjointe 
p{f,w,Tie,9f,cx;g) et la loi a posteriori p{f ,w, 11^,9 f,cx\g). Cependant l'expression de 
cette loi 

P{f,w,'^e,0f,cx\g) oc p{g\w, f p{f\9f) p{w\a) 

p{ee)p{9f)p{a) 

n'est pas séparable en ses composantes. L'approche variationnelle consiste donc à l'appro- 
cher par une loi séparable 

p{f, w, de, 9f, a\g) ~ q{f)q{w) [] g(^ejg(%) H (^4) 

i 3 

et avec les choix des lois a priori conjuguées en appliquant la procédure décrite plus haut, 
on obtient 

q{f) = ■A^(M/,S/) avec 
S/ = < W^BW >^+<9f> D}Df]-\ 
Hf = S/** <W>' Bg, (45) 

q{w) = AfiHy,, S^) avec 

= [^^ <F^BF>^ + A]~\ 
fj,^ = ^^^'<F>'Bg (46) 

g(6'ej = Qi(yei,Pei) aVCC 
«ei «60 + -^/2, 

(47) 

g(6'/) = Ç{af,Pf) avec 
a/ = afo + N/2, 

Pf = /?yo + 1/2 Tr{£)}£)^ <//*>}, (48) 
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qi^wj) = G{o,j,bj) avec 
Qj = a^o + 1/2, 

bj = /3^o + 1/2 < > . (49) 

où 

A = diag [< a^j- >, j = 1, • • • , A^] et 

B-diag[</3ei >,i-l,-- - , ^] ■ 
On a ainsi l'expression des différentes composante de la loi scparable approchante. On 
peut en déduire facilement les moyennes de ces lois, car ces lois sont, soit des gaussiennes, 
soit des lois gamma. 

<w>= ix^, < w] >= [ixj] + [S^]^-j- 

<f>^fif, <//*>=/x^/x* +S/ (50) 

< 9ei >= aei/Pei < Of >= Oif/Pf, 

et 

< ee* > = oo* 2g\< F > ^ < w >Y 
<w><w>'^ F*]** 

Pour le calcul des termes < W*W > et < F*F > qui interviennent dans les expres- 
sions de S/, et [Fww^F^] on peut utiliser le fait que F et W sont des matrices 
bloc-Toeplitz avec des blocs Toeplitz (TBT) , on peut les approcher par des matrices bloc- 
circulantes avec des blocs circulantes (CBC) et les inverser en utilisant la TFD. Notons 
aussi que Hf et peuvent être obtenu par optimisation de 

Jil^f) = [g- ^WiXfYB[g - <0f> \\Dfl^f\? (52) 

et de 

J{^J^w) = [9- ^FnjB[g - ^FuJ + \\AiJ,J\ (53) 



5 Restauration avec modélisation Gauss-Markov-Potts 

Le cas d'une modélisation gaussienne reste assez restrictif pour la modélisation des 
images. Des modélisation par des champs de Markov composites (intensités-contours ou 
intensités-régions) sont mieux adaptées. Dans ce travail, nous examinons ce dernier. L'idée 
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de base est de classer les pixels de l'imagcs / = {/(r),r G TZ} en K classes par l'in- 
termédiaire d'une variable discrète z{r) e {1, • " " )-^}- L'image z — {z{r),r G TZ} 
représente ainsi la segmentation de l'image /. Les pixels fk — {f{r),r G TZ^} avec 
TZk — {r : z{r) — k} ont des propriétés communes, par exemple la même moyenne /ik et 
la même variancc (homogénéité au sens probabiliste) . Ces pixels se trouvent dans un 
nombre fini de régions compactes et disjointes TZki telles que : UiTZki = TZk et UkTZk = TZ. 
On suppose aussi que fk et fi^k ^ l sont indépendants. 

A chaque région est associée un contour. Si on représente les contours de l'images 
par une variable binaire c(r), on a c(r) = à l'intérieure d'une région et c(r) = 1 aux 
frontières de ces régions. On note aussi que c(r) s'obtient à partir de z{r) d'une manière 
déterministe (voir Fig 1). 

Avec cette introduction, nous pouvons définir 

p{f{r)\z{r) = k, rrik, Vk) = N{mk, Vk) (54) 
ce qui suggère un modèle de mélange de gaussienne pour les pixels de l'image 

= ^o.kM{mk,Vk) avec = P{z{r) = k) (55) 

k 

Une première modélisation simple est de supposer que les étiquettes z{r) sont a priori in- 
dépendants : 

r 

Nous appelons ce modèle. Mélange séparable de gaussiennes (MSG). 

Maintenant, pour prendre en compte la structure spatiale de ces pixels, nous devons 
introduire, d'une manière ou autre, une dépendance spatiale entre ces pixels. La modéli- 
sation markovienne est justement l'outil approprié. 

Cette dépendance spatiale peut être fait de trois manières. Soit utiliser un modèle 
markovien pour z{r) et un modèle indépendant pour f{r)\z{r), soit un modèle markovien 
pour f{r)\z{r) et un modèle indépendant pour z{r), soit un modèle markovien pour 
f{r)\z{r) et un modèle markovien aussi pour z{r). Nous avons examiné ces cas avec des 
modèles de Gauss-Markov pour f{r)\z{r) et le modèle de Potts pour z{r). Ce dernier 
peut s'écrire sous deux formes : 



p{z{r)\z{r'),r' G V(r)) oc exp 



7 Y. Kz{r)-z{r')) 



(57) 
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/(r) z{r) c{r) 

FiG. 1 - Modèle de mélange et champs de Markov caché : image des intensités ou niveau 
de gris f{r), image z{r) de segmentation ou classification, image binaire c(r) des contours. 



et 



p{z) oc exp 

Ces différents cas peuvent alors se résumer par : 



7E E à{^ir) - z{r')) 

reTlr'eVir) 



(5^ 



Modèle de mélange séparables de gaussiennes (MSG) : 

C'est le modèle le plus simple oii aucune structure spatiale n'est pris en compte a priori . 
Les relations qui donnent la loi a priori conjointe de z) = p{f\z) p{z) sont : 



p{f{r)\z{r) = k) 
p{f\z) 



avec rrizir) = m^, Vr G TZk et Vz{r) = f^, Vr G TZk, et 

ak, Vr G 7^ 



(59) 




a 



(60) 



avec rik = J^ren ^(-2(^) — k) le nombre de pixels dans la classe k et = n le nombre 

total des pixels de l'image. Les paramètres de ce modèle sont 

= {(«fc, mk,Vk),k = 1, - ■ ■ , K}. 

Rappelons aussi que ^^^a/s = 1. Lors d'une estimation bayésienne non supervisée, il faut 
aussi attribuer des lois a priori à ces paramètres, les lois conjuguées correspondantes sont 
des gaussiennes pour m^, des Inverse Gammas pour Vk et la loi Dirichlet pour et = [ak, k = 
,K] : 
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p{mk\mo,Vo) = J\f{mo,Vo), Mk 

p(wfc|ao,/3o) = IG{ao,bo), \/k (61) 

p(o:|ao) = î^(aolK) 

où mo, f 0, Oo; ^0 et «o sont fixés pour un problème donnée. On les choisis de telle sorte que 

ces lois soient les moins informatives (par exemple, mo = 0, fo = 10, oq = 1, 6o = 100 et 

«0 = l/K). 

Modèle de mélange séparable de Gauss-Markov (MSGM) : 

Ici, la structure spatiale est pris en compte au travers d'un modèle markovien sur / : 



p(/(r)|^(r), /(rO, z(r'), r' G V(r)) = U{m,{r)Mr)) 



(62) 



avec 



|V(r)| ^r'eV(r) ^z(r') 

fi^(r') si z{r') ^ z{r) (63) 
f{r') si z{r') = z{r) 
Par contre le champs z est supposé séparable comme dans le premier cas (160]) . 

On remarque que f\z est un champ de Gauss-Markov non homogène car la moyenne 
mz{r) et la variance Vz{r) varient en fonction de r. Tous les pixels fk = {/(r),r G TZk} 
se trouvant dans une région k forment alors un vecteur gaussien de moyen rriklk et de 
matrice de covariance : 

p{fk) = {rriklk, VkTik) (64) 

où ffcSfc est une matrice de covariance de dimension x et 1^ = 1 est un vecteur de 
taille Uk remplis de 1. Par ailleurs, on a aussi : 



m 



zir') 



P{f\z) OC Y\_-^{rnklk,Vk'Sk) 



(65) 



ce qui signifie que les pixels appartenant aux différentes régions sont indépendantes. 



Modèle de Gauss-Potts (MGP) : 

Dans ce modèle, la structure spatiale est pris en compte au travers du modèle markovien 
de Potts (l57j) ou (|58|1 . La loi de p{f\z) est la même que dans le premier cas. En résumé, 
ici on a : 

' Pim = Uren^{rn{r),v{r)) 



p{z) 



oc exp 



7 J2ren Er'ev(r) ^i^r) - z{r')) 



19 



zCr-> 

MSG 



► z(r) 

MSGM 



MGP 




MGMP 



FiG. 2 - Quatre modèles de mélange avec champs cachés. Dans le modèle MSG toutes les 
variables sont séparables en r. Dans MSGM, /(r) dépends de z{r) et de ses voisins f{r'). 
Dans MGP, z{r) dépends de ses voisin z{r'). Dans MGMP, z{r) dépends de ses voisin 
z{r') et /(r) dépends de z{r), z{r') et de ses voisins f{r'). 

avec m(r) = mk,'ir G IZk et v{r) = ffc, Vr G 7?.^. 

Ici, donc, les pixels de l'images sont indépendantes conditionnellement à la connais- 
sance des variables cachées, et donc, contrairement au cas précédent (!64|) . ici on a p{fk) = 
A/'(mfclfc, ffc/fc) où ffc/fc est une matrice de covariance identité de dimension x n^, ce 
qui permet d'écrire : 



P{fk) = Uiniklk^Vklk) oc exp 



-E 



f{r) - mk) 

Vk 



(67) 



et 



p{f\z) oc JjA/'(mfclfc, Ifc) 



oc exp 



(68) 



k r€Tlk 

ce qui signifie que, comme le cas précédent, les pixels appartenant aux différentes régions 
sont indépendantes. 

Les paramètres de ce modèle sont Of = {(«fc, f^k, Vk), k = 1, - ■ ■ , K} et 7. Malheureu- 
sement, il n'y a pas de loi conjuguée pour ce paramètre et donc son estimation devient 
plus difficile. Dans ce travail, nous fixons ce paramètre a priori . 

Modèle de Gauss-Markov-Potts (MGMP) : 

Il s'agit là de la composition des deux derniers modèles. Ici, nous résumons ses relations 
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d'une manière légèrement différente en utilisant : 



Crir') = 1 - ô{z{r') - z{r)) 



1 si z{r') 7^ z{r) 
si z{r') = z{r) 



(69) 



ce qui donne : 



pifirMr), fir'), ^(r'), r' G V(r)) = M{m,{r),v,{r)) 



(70) 



avec 



m^(r) 



1 



J2 [(l-c,(r'))^.(o + c.(r')/(r')] 



(71) 



|V(r)| 



r'eV(r) 



et 



p(2;)ocexp XI (l-c,.(r')) 



(72) 



re7^r'GV(r) 

En tous cas, quelque soit le modèle choisi parmi ces différents modèles, l'objectif est 
d'estimer f,zetO en utilisant la loi a posteriori jointe : 



Bien que nous connaissons les expressions de tous les composants de numérateur de la 
fraction à la droite de cette relation, le calcul du dénominateur n'est que rarement possible 
d'une manière analytique. On cherche alors à approcher p{f,z,0\g) par un produit des 
lois plus simples à manipuler. Un premier choix est de l'approcher par une loi séparable 
q{f,z,6\g) = qi{f) 92(2) Qsi^)- Un deuxième choix est qi{f\z) q2{z) q-s^O) qui permet 
de garder des liens forts qui existent entre / et 2; et de ne relaxer que des liens faibles 
entre 6 et ces deux derniers. 

Le choix des familles appropriées pour qi{f\z) et q2{z) et qsi^O) pour chacun de ces 
modèles, qui est en lien avec les formes des lois a priori dans chacun de ces cas, et 
les expressions de la mise à jour de ces différentes lois au cours des itérations nécessite 
beaucoup d'espace. Ici, nous allons juste fournir le principe et un résumé de ces relations. 

Le choix de qsiO), ou plus exactement, des familles des lois pour chacune des compo- 
sentes de 6, qui sont la même pour tous les cas, sont des lois conjuguées de flïïTj) . Par la 
suite, nous détaillons les choix de q{f\z) et q{z) dans les différent cas. 




(73) 
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(74) 



Cas 1 (MSG) : 

Dans le cas du premier modèle on a 

PiM = UMfirMr)) 

OÙ p{f {r)\z{r) = k) = J\f{'mk,Vk) et p{z{r) = k) = ak, ce qui naturellement nous conduit 
à choisir 

' qif\z) = UrQifirMr)) 

avec q{f{r)\z{r) = k) ^J\f{mk,Vk) et q{z{r) = /c) = à^. 
Au cours des itérations, Vk et àk seront mise à jour. 

Cas 2 (MSGM) : 

Dans le deuxième cas, le modèle de mélange de Gauss-Markov, notant que 

P{f\z) = l[m{r)\m,ir),v,{r)) 



ou 



avec 



1 V- * 



(75) 



(76) 



(77) 



(78) 



m:(,,) = 5{z{r') - z{r))f{r') + (1 - 5{z{r') - z{r))^x,^r') 
on trouve naturellement l'approximation suivante : 

m:(,,) = 5{z{r') - z{r))f{r') + (1 - 5{z{r') - ^(r))m,M (79) 

où f{r') est l'espérance de f{r') calculée à l'étape précédente. Il s'agit de l'approximation 
dite en champs moyens pour le champ de Markov / : 



ç(/l^)=n.g(/(r-)|^(r-)) = Y{r^{î{r)\m,{r)Mr)) 
avec 

^z{f) = |v(r) Xlr'GV(r) 

5{z{r) - z{r'))f{r) + (1 - S{z{r) - ;^(r')))/x.w 
Cas 3 (MGP) : 

Dans le troisième cas, le modèle de Gauss-Potts on a 

Pim = UrPifir)\z{r)) 
P{z) = Y{rP{^{r)\z{r')y eV{r)) 



(80) 



(81) 
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ce qui naturellement nous conduit à choisir 

= Ur<lifir)\zir)) 




(82) 

Ur<li4rMr'),r' eV{r)) 

où nous avons choisi une approximation en champs moyens pour le champs de Potts. 
Cette approximation supplémentaire qui consiste à remplacer les valeurs de z{r') par 
leurs moyenne z{r') est une approximation courante pour un modèle de Potts. 

Cas 4 (MGMP) : 

Dans le quatrième cas,le Modèle de Gauss-Markov-Potts on a 
avec 

ô{z{r) - z{r')f{r) + (1 - ô{z{r) - z{r'))fi,^r) 

Dans tous les cas, puisque les lois a priori pour les différent composants (\e d — 
{{^ej}) {"^fe}) {vk\-i {Pk}} sont des lois conjuguées et séparables, on choisi 

m = n ^(^^) (84) 

i 

Il reste un paramètre que nous garderons fixe. C'est le paramètre du modèle de Potts. 
En effet, n'ayant pas une expression analytique pour la fonction de répartition du modèle 
de Potts, il n'existe pas une loi conjuguée pour ce paramètre. Dans un premier temps 
donc, nous gardons ce paramètre comme un paramètre de réglage non supervisé. 

Les deux tableaux qui suivent résument les lois a priori et les choix des lois séparables 
pour les quatre cas proposés. 

Ici, nous ne développons pas plus les expressions de la mise à jour de ces différentes 
loi au cours des itérations, mais ce qu'il faut savoir est que toutes ces lois étant de formes 
paramétriques connues (gaussienne, gamma ou inverse gamma, Wishart ou inverse Wi- 
shart et Dirichlet), nous avons des expressions analytiques pour les espérances de ces lois. 
Les détails de ces relations accompagnés des résultats de simulation seront communiqués 
dans un autre papier dans un future proche. 
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Modèle 


MSG 


MSGM 


MGP 


MGMP 


P{f\z) 


WM{m^{r),v^{r)) 

r 


Gauss-Markov 


YlAf{m,{r),v,{r)) 

r 


Gauss-Markov 


p{z) 


r k 


r k 


rOXXS 


r OttS 


p{mk) 


J\f{mo,vo) 


J\f{mo, vo) 


J\f{mo,vo) 


J\f{mo,vo) 






XÇ{ao,bo) 




1Q{ao, bo) 


p{ct) 


©(«o, ■ ■ ■ ,"o) 




V^ao, ■■■ ,«o) 


V{ao, ■■■ ,ao) 


P{0e) 


S{aeo:beo) 


S{aeo:beo) 


G{aeo:beo) 


G{aeo,beo) 




r 


r 


Y[Afim^ir),v,{r)) 

r 


r 


q{z) 


r k 


r k 


Ylq{z{r)\z{r'),r'eV{r)) 

r 


r 


q{mk) 


M{mk,Vk) 


M{mk,Vk) 








W{àk,bk) 


W{àk,bk) 


IG{àk,bk) 


W{àk,bk) 


q{a) 


T>{âi, ■ ■ • ,âx) 


T>{âi, ■ ■ • ^Ûk) 


T>{âi, ■ ■ ■ ,âx) 


T>{ài, ■ • ■ ,âx) 


q{Oe) 


Giàeo,beo) 


Giàeo,beo) 




^(âeo,^€o) 



Tab. 1 - Les lois a priori et les lois séparables approchantes pour les quatre modèles 
proposés. 
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6 Résultats et discussions 



Nous avons utilisé cette approche dans plusieur domaines de problèmes inverses : i) 
restauration d'image [53l|42], ii) reconstruction d'image en tomographie X [541 155], et iii) 
séparation de sources et segmentation des images hyperspectrales [3D] . Nous sommes aussi 
en train de mettre en oeuvre ces algorithmes pour le cas 3D de la tomographie X ainsi 
qu'au cas de l'imagerie microondes. 

Dans tous ces applications, nous avons implémenté à la fois des algorithmes MCMC 
(échantillonneurs de Gibbs) et l'approche variât ionnelle bayésienne (VB) développée dans 
cet article. Nous avons ainsi pu comparé les résultats obtenus par ces deux méthodes. 
D'une manière générale, les principales conclusions de ces expérimentations et compa- 
raisons sont : i) la qualité des résultats obtenus sont pratiquement similaires, et ii) le 
principal avantage de l'approche VB est dans le gain en temps de calcul qui dépend bien 
sure du contexte et de l'applications. 

Un autre avantage important que l'on peut mentionner est l'existance d'un critère 
(l'énergie libre, équation [TU]) que l'on peut utiliser : i) comme critère d'arrêt pour l'algo- 
rithme et ii) comme un critère de choix de modèle. En effet, comme nous l'avons vu, grâce 
à la relation in minimiser KL(g : p) est équivalent à maximiser J-'{q) et sa valeur optimale 
est un bon indicateur pour \np{g\A4). Ainsi, les valeurs relatives de J-'{q) au cours des 
itérations de l'algorithme peuvent être utilisées comme un critère d'arrêt et sa valeur op- 
timale atteint pour un modèle Aii peut être comparée à sa valeur optimale atteint pour 
un autre modèle Ai2 comme un critère de préférence entre les deux modèles. 

7 Conclusion 

L'approche variationnelle de l'approximation d'une loi par des lois séparables est ap- 
pliquée au cas de l'estimation non supervisée des inconnues et des hyperparamètres dans 
des problèmes inverses de restauration d'image (déconvolution simple ou aveugle) avec des 
modélisations a priori gaussiennes, gaussiennes généralisées, mélange de gaussiennes indé- 
pendantes ou mélange de Gauss-Markov avec champs cachés des étiquettes indépendantes 
ou markovien (Potts). 

Dans ce papier, nous avons décrit d'abord le principe de l'approximation d'une loi et 
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ensuite proposé des approximations appropriées pour des lois a posteriori conjointes que 
l'on trouve dans le cas des problèmes inverses en général et en restauration d'images en 
particulier. Cependant, la mise en oeuvre effective de ces méthodes est en cours et les ré- 
sultats de simulation et évaluation des performances de ces méthodes seront communiqués 
dans un future proche. 

Les détails de ces relations accompagnés des résultats de simulation qui sont disponible 
sous forme d'un rapport [53l |42] seront communiqués dans un autre papier en préparation 
dans un futur proche. 
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L'author souhaite remercier Jean- François Berchet et Guy Demoment pour les discus- 
sions fructueuses et la relecture de la première version de ce papier, H. Ayasso et S. Fekih 
pour l'implémentation de ces algorithmes en restauration d'image et en tomographie X, 
N. Bali pour l'implémentation de ces algorithmes en séparation de sources et en imagerie 
hyperspectrale. 

Références 

[1] J. Rustagi, Variational methods in statistics. New York : Académie Press, 1976. 

[2] Z. Ghahramani and M. Jordan, "Factorial Hidden Markov Models," Machine Learning, 
no. 29, pp. 245-273, 1997. 

[3] W. Penny and S. Roberts, "Bayesian neural networks for classification : how useful is the 
évidence framework?," Neural Networks, vol. 12, pp. 877-892, 1998. 

[4] S. Roberts, D. Husmeier, W. Penny, and I. Rezek, "Bayesian approaches to gaussian mix- 
ture modelling," IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 20, 
no. 11, pp. 1133-1142, 1998. 

[5] H. Attias, "Independent factor analysis," Neural Computation, vol. 11, no. 4, pp. 803-851, 
1999. 

[6] M. Jordan, Z. Ghahramani, T. Jaakkola, , and L. Saul, "An introduction to variational 
methods for graphical models," Machine Learning, vol. 37, pp. 183-233, 2006. 



26 



[7] W. Penny and S. Roberts, "Dynamic models for nonstationary signal segmentation," Com- 
puters and Biomédical Research, vol. 32, no. 6, pp. 483-502, 1999. 

[8] H. Attias, "A variational Bayesian framework for graphical models.," in Advances in Neural 
Information Processing Systems (S. SoUa, T. K. Leen, and K. L. Muller, eds.), vol. 12, 
pp. 209-215, MIT Press, 2000. 

[9] T. Jaakkola, "Tutorial on variational approximation mcthods," in Advanced mean field me- 
thods : theory and practice. (M. Opper and D. Saad, eds.), (Cambridge, Massachusetts), 
pp. 129-159, 2000. 

[10] J. Miskin, Ensemble Learning for Independent Component Analysis. Thèse de doctorat, 
Cambridge, 2000, http : //www. inference.phy.cain.ac.uk/jwinl003/. 

[11] J. W. Miskin and D. J. C. MacKay, "Application of ensemble learning to infra-red imaging," 
in Proceedings of the Second International Workshop on Independent Component Analysis 
and Blind Signal Séparation, pp. 399-404, 2000. 

[12] J. W. Miskin and D. J. C. MacKay, "Ensemble learning for blind source séparation," in 
ICA : Principles and Practice (S. Roberts and R. Everson, eds.), Cambridge : Cambridge 
University Press, 2001. 

[13] W. Penny and S. Roberts, "Bayesian multivariate autoregresive models with structured 
priors," lEE Proceedings on Vision, Image and Signal Processing, vol. 149, no. 1, pp. 33-41, 
2002. 

[14] S. Roberts and W. Penny, "Variational bayes for generahsed autoregressive models," IEEE 

Transactions on Signal Processing, vol. 50, no. 9, pp. 2245-2257, 2002. 

[15] M. Cassidy and W. Penny, "Bayesian nonstationary autogregressive models for biomédical 
signal analysis," IEEE Transactions on Biomédical Engineering, vol. 49, no. 10, pp. 1142- 
1152, 2002. 

[16] W. Penny and K. Priston, "Mixtures of gênerai linear models for functional neuroimaging," 
IEEE Transactions on Médical Imaging, vol. 22, no. 4, pp. 504-514, 2003. 

[17] W. Penny, S. Kiebel, and K. Priston, "Variational Bayesian inference for fmri time séries," 
Neurolmage, vol. 19, no. 3, pp. 727-741, 2003. 

[18] R. A. Choudrey and S. Roberts, "Variational bayesian mixture of independent component 

analyscrs for finding sclf-similar arcas in images," in Proc. 4th International Symposium on 
Independent Component Analysis and Blind Signal Séparation (ICA2003), (Nara, Japan), 
pp. 107-112, April 2003. 



[19] R. Choudrey and S. Robcrts, "Variational Mixture of Bayesian Independent Component 
Analysers," Neural Computation, vol. 15, no. 1, 2003. 

[20] W. Penny, R. Everson, and S. Roberts, "Hidden markov independent component analysis," 
in Advances in Independent Component Analysis (M. Giroliami, éd.), Springer, 2000. 

[21] D. M. Blei and M. I. Jordan, "Variational methods for the dirichlet process," in ICML, 2004. 

[22] N. Nasios and A. Hors, "A variational approach for Bayesian blind image deconvolution," 
IEEE Transactions on Signal Processing, vol. 52, no. 8, pp. 2222-2233, 2004. 

[23] C. Archambeau, T. Butz, V. Popovici, M. Verleysen, and J. Thiran, "Supervised nonparame- 
tric information theoretic classification," in Proceedings of the 17th International Conférence 
on Pattern Récognition, vol. 3, pp. 414-417, IEEE, August 2004. 

[24] M. Woolrich, E. Timothy, F. Christian, and S. Smith, "Mixture models with adaptive spatial 
regularization for segmentation with an application to fmri data," IEEE Trans. on Médical 
Imaging, vol. 24, pp. 1-11, January 2005. 

[25] D. M. Blei and M. I. Jordan, "Variational methods for the dirichlet process," Journal Version 
(TODO : which ?), 2006. 

[26] M. Beal and Z. Ghahramani, "Variational Bayesian learning of directed graphical models 
with hidden variables," Bayesian Statistics, vol. 1, pp. 793-832, 2006. 

[27] H. Kim and Z. Ghahramani, "Bayesian gaussian process classification with the em-ep algo- 
rithm," IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 28, no. 12, 
pp. 1948-1959, 2006. 

[28] B. Wang and D. Titterington, "Convergence properties of a gênerai algorithm for calculating 
variational bayesian estimâtes for a normal mixture model," Bayesian Analysis, vol. 1, 
pp. 625-650, 2006. 

[29] K. Watanabe and S. Watanabe, "Stochastic complexities of gaussian mixtures in variational 
bayesian approximation," Journal of Machine Learning Research, pp. 625-644, April 2006. 

[30] N. Nasios and A. Bors, "Variational learning for gaussian mixture models," IEEE Transac- 
tions on Systems, Man and Cybemetics, Part B, vol. 36, no. 4, pp. 849-862, 2006. 

[31] K. Friston, J. Mattout, N. Trujillo-Barreto, J. Ashburner, and W. Penny, "Variational free 
energy and the laplace approximation," Neuroimage, no. 2006.08.035, 2006. Available On- 
line. 



28 



[32] W. Penny, S. Kicbcl, and K. Friston, "Variational bayes," in Statistical Parametric Mapping : 
The analysis of functional brain images (K. Priston, J. Ashburner, S. Kiebel, T. Nichols, 
and W. Penny, eds.), Elsevier, London, 2006. 

[33] Z. Ghahramani, T. Griffiths, and P. SoUich, "Bayesian nonparametric latent feature models," 
Bayesian Statistics, vol. 8, 2007. 

[34] C. McGrory and D. Titterington, "Variational approximations in Bayesian model sélection 
for finite mixture distributions," Computational Statistics and Data Analysis, vol. 51, no. 11, 
pp. 5352-5367, 2007. 

[35] F. Forbes and F. Gersende, "Combining monte carlo and mean-field-like methods for infe- 
rence in hidden markov random fields," IEEE Trans. on Image Processing, vol. 16, pp. 824- 
835, March 2007. 

[36] M. Ichir and A. Mohammad-Djafari, "A mean field approximation approach to blind source 
séparation with Ip priors," in Eusipco 2005, Antalya, Turkey, September 2005, Eusipco 2005, 
Antalya, Turkey, September 2005, September 2005. 

[37] M. Ichir and A. Mohammad-Djafari, "Hidden markov models for wavelet image séparation 
and denoising," in Proc. IEEE International Conférence on Acoustics, Speech, and Signal 
Processing (ICASSP '05), vol. 5, pp. v/225-v/228, 18-23 March 2005. 

[38] H. Snoussi and A. Mohammad-Djafari, "Estimation of structured gaussian mixtures : The 
inverse em algorithm," IEEE Trans. on Signal Processing, vol. 55, pp. 3185-3191, July 2007. 

[39] N. Bali and A. Mohammad-Djafari, "Mean Field Approximation for BSS of images with 
compound hierarchical Gauss-Markov-Potts model," in MaxEnt05,San José CA, US, Ame- 
rican Institute of Physics (AIP), August 2005. 

[40] N. Bali and A. Mohammad-Djafari, "Bayesian approach with hidden markov modeling 
and mean field approximation for hyperspectral data analysis," IEEE Trans. on Image 
Processing, vol. 17, pp. 217-225, Feb. 2008. 

[41] A. Mohammad-Djafari, "Approche variationnelle pour le calcul baycsien dans les problèmes 
inverses en imagerie," in GRETSI, Troyes, France, September 2007. 

[42] H. Ayasso and A. Mohammad-Djafari, "Variational bayes with gauss-markov-potts prior 
models for joint image restoration and segmentation," Visapp Proceedings, (Funchal, Ma- 
daira, Portugal), Int. Conf. on Computer Vision and Applications, 2008. 



29 



[43] A. Mohammad-Djafari, "Gauss-markov-potts priors for images in computer tomography 
resulting to joint optimal reconstruction and segmentation," International Journal of To- 
mography & Statistics, vol. 11, no. W09, pp. 76-92, 2008. 

[44] R. Molina, A. K. Katsaggelos, and J. Mateos, "Bayesian and regularization methods for 
hyperparameter estimation in image restoration," IEEE Transactions on Image Processing, 
vol. 8, pp. 231-246, February 1999. 

[45] A. Likas and N. Galatsanos, "A variational approach for Bayesian blind image deconvolu- 
tion," IEEE Trans. on Signal Processing, vol. 52, pp. 2222-2233, August 2004. 

[46] K. Blekas, A. Likas, N. Galatsanos, and I. Lagaris, "A spatially-constrained mixture model 
for image segmentation," IEEE Trans. on Signal Processing, vol. 16, pp. 494-498, March 
2005. 

[47] A. Mohammad-Djafari and L. Robillard, "Hierarchical Markovian models for 3D computed 
tomography in non destructive testing applications," in EUSIPCO 2006, EUSIPCO 2006, 
September 4-8, Florence, Italy, September 2006. 

[48] M. Patriksson, Nonlinear programming and variational inequality problems. A unified ap- 
proach. Applied Optimization, Dordrecht, The Netherlands : Kluwer Académie Publishers, 
May 1999. 

[49] R. Choudrey, W. Penny, and S. Roberts, "An ensemble learning approach to independent 
component analysis," in IEEE Workshop on Neural Networks for Signal Processing, Sydney 
Australia, 2000. 

[50] W. Penny and S. Roberts, "Variational bayes for non-gaussian autoregressive models," in 
IEEE Workshop on Neural Networks for Signal Processing, Sydney Australia, 2000. 

[51] R. Choudrey and S. Roberts, "Variational Bayesian Mixture of Independent Component 
Analysers for Finding Self-Similar Areas in Images," in ICA, NARA, JAPAN, April 2003. 

[52] B. R. Hunt, "A matrix thcory proof of the discrète convolution theorem," IEEE Transactions 
on Automatic and Control, vol. AC-19, pp. 285-288, 1971. 

[53] H. Ayasso and A. Mohammad-Djafari, "Approche bayésienne variationnelle pour les pro- 
blèmes inverses, application en tomographie microonde," tech. rep., Rapport de stage Master 
ATS, Univ Paris Sud, L2S, SUPELEC, 2007. 

[54] H. Ayasso, S. Fekih-Salcm, and A. Mohammad-Djafari, "Variational bayes approach for to- 
mographie reconstruction," in Proceedings of the 28th International Workshop on Bayesian 



30 



Inference and Maximum Entropy Methods in Science and Engineering, MaxEnt, vol. 1073, 
(Boraceia, Sao Paulo (Brazil)), pp. 243-251, AIP, November 2008. 

[55] H. Ayasso and A. Mohammad-Djafari, "Joint image restoration and segmentation using 
gauss-markov-potts prior models and variational bayesian computation," Submitted to IEEE 
Image Processing, Febuary, 2009. 



31 



